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CHAPTER-1 

INTRODUCTION 


1.1 Definition 

A thermosyphon is a circulating fluid system 
whose motion is caused by density differences in a body 

force field which result from heat transfer .This definition 
given by Davies and MorrisCID is so broad as to include 

all natural convection processes plus others. The system to 

which the name thermosyphon is applied in formal studies 
are in fact systems which has the intrinsic function of 
removing heat from a prescribed source and transporting 
heat and mass over a specific path ( frequent ly a 
reci r culat ing flow)and rejecting the heat and mass to a 

prescribed sink. Thus the path of the recirculating flow 
which transports the thermal energy is or can be totally 

prescribed. Thermosyphon flows are intrinsically driven by 
thermal buoyancy forces either locally or in an overall 

sense. A simple loop flow may be the result of local 

buoyancy force alone, but a mult ibranched flow circuit can 
incorporate sections in which the flow direction is 

contrary to local buoyancy force resulting from pressures 
created by the overall system buoyancy forces. Based on 

these factors.a thermosyphon is defined as a prescribed 

circulating fluid system driven by thermal buoyancy 
forces . 

1.2 Classification 

The thermosyphons can be classified into 
the following categories according to 

a) The nature of boundaries! 
i)Open thermosyphon 

ii)Closed thermosyphon 

b) The regime of heat transfer: 
i)Purely natural convection 

ii)Mixed natural and forced convection 

c) The number or types of phases present: 
i)Single phase 

ii)Two phase 



d)The nature of the body forces 
i ) Gravitational 
i i IRotational 
1 .3 Application 

The most common industrial thermosyphon 

applications include the following: 

i)Gas turbine blade cooling 
ii lElectrical machine rotor cooling 

iii )Transf ormer cooling 

iv)Nucfear reactor decay heat removal 

v)Heat exchanger fins 

vi)Cryogenic cool-down apparatus 

vii)Steam tubes for bakers' ovens 
viiilCooling of internal combustion engines 

ix)Solar water heater 
xJConvection in the earth's mantle 
xi JTemperature distribution in earth drillings in steam 
power fields 

xii)Preservations for the conservation of permaforst under 


buildings in the Canadian Northland 
xi i i IMaintainance of ice-free navigation buoys 




A 


variety 

of 

thermosyphon 

characteristics 

are 

responsible for 

the applications found 

till now and 

can 

lead 

to 

numerous 

further 

applications. 

i)A thermosyphon 

can 

act 

as a 

thermal 

conductor with 

either a 

smal 1 

or 

a 

large thermal impedance. 

ii)It can be 

used 

as 

a 

thermal 

diode 

or rectifier , or! 


even as a thermal triode, permitting a variation in heat 1 

flow based on small changes in temperature. 

1.4 Decay Heat Removal in Nuclear Power Plant 

In a nuclear reactor core , enormous! 
amount of heat is generated in the fuel element due to; 
the fission of the fissionable material present in it. This! 

heat is taken out of the reactor core by the coolant by! 

forced convection by means of primary pumps. The primary; 

! 

coolant may exchange heat to the working fluid of the; 
turbine-condensfcr combination in the heat exchanger.lt may j 
be the working fluid itself as in the direct-cycle BWR.Ini 



normal operating conditions, the coolant flow rate and the 
coolant inlet temperature are so designed as to take out 

the heat from the fuel element in such a manner that the 

fuel element is not damaged due to excessive temperature. 

However, a nuclear reactor is not 

only designed to work under normal conditions , but also it 
is to be seen that under accident condi t ions , the safety 

limits are not exceeded. As the nuclear industry came into 
recognition after a much horrified and fearful occurence 
during the Second World war, the general faith of the 
ordinary public is very low with it. Again to avoid the 

dangers involved with the release of radio— activity, much 

care is taken in the design, construction, planning and 

operation of a nuclear power plant both under normal and 

accident conditions. 

In nuclear reactors , loss of flow 

accidents (LOFA) is one of the major accidents , caused by a 
flow obstruction in the core or because of a decrease of 

pressure drop due to a pump failure. The pump may fail 

due to mechanical seizure or electrical supply failure. In 

the former case, the flow may be assumed to reduce to 

zero instantaneously whereas in the latter case, the flow 
deteriorates rather gradually due to the high inertia of 

the flywheels. 

During the normal operations , fuel 

temperatures are controlled by the forced circulation of 
the coolant. However , during certain events postulated to 

occur during the life-time of a reactor, power to the 

primary pump, which maintains the forced flow, may be lost. As 
the forced flow tends to decrease, the plant protective 
system will, with extremely high rel iabi 1 ity, shut down the 
plant .However , due to the generation of the radioactive 

material within the fuel and the structural components of 

the reactor during power operation , residual decay heat will 
be there even though the primary fission reactions have 
ceased. Theref ore, there must be some continued circulation of 

coolant , though at much reduced rate , fol lowing this event to 

prevent overheating of the reactor. 



Although the design of the plant 
will include back-up pumping systems such as redundant 
power supplies, the ultimate fall-back heat transport 
mechanism is natural circulation cooling. Provided the 

reactor has been shut down fast enough to prevent damage 

during the flow coast-down, a central question then becomes 
whether an equilibrium state can be established f or the 

natural convection removal of decay heat without the 
occurence of excessive temperatures. 

The natural circulation decay 

heat removal in PHWR is discussed now. The PHWR has several 
modes of decay heat removal. The boilers and the primary 
pumps provide the normal mode. The shutdown cooling system 
has independent heat exchangers and pumps and can remove 
decay heat at pressures upto nominal primary system 
pressures .Thermosyphoning is the mode of heat transport to 

the boilers in the absence of any forced circulation of 
the primary coolant. 

Thus thermosyphoning is defined as the 
natural convective flow of the primary coolant over the 

boilers.lt is the predicted mode of heat transport in 
many postulated scenarios for PHWR safety analysis.The 
scenarios encompass a wide range of primary coolant 

inventories and secondary side conditions. 

A typical figure-of-eight loop of a 
PHWR primary heat transport system is shown in Figure-1. It 
comprises two passes of the coolant past the fuel. Thus 

the coolant sees Corel , boi le rl , pumpl , core2, boi lerE, pump2, Corel 
etc.. Each core pass includes a large number of horizontal 
fuel channels each of which is connected by 
an outlet feeder to large diameter inlet 
headers .Large diameter pipings are connecting 
to the boilers and pumps. The tops of boiler 
located typically 11m. above the headers and 
are located 7m. above the channels. The total 


an inlet and 
and outlet 
the headers 
U-tubes are 
the headers 
height is 


18m. and the density differnces in this height provide 
the buoyancy force for thermosyphoning. 

PHWR primary pumps are provided 



of electrical 


with high inertia flywheels which on loss 
supply , provide a forced flow of coolant whilst power is 

reduced to decay level .Subcool ed ( single— phase ) thermosyphoning 
is predicted thereafter for this and the more probable 

accident scenarios. 

Boi 1 ing ( two— phase ) thermosyphoning is 

predicted either under reduced primary coolant inventory 
conditions or if the coolant make-up system is inoperative 
during boiler depressur isat ion . In the latter case,a 
• r .-w'.i' ! ri;er compensates , in part, for shrinkage of the primary 
coolant as it cools down. In any event , two-phase 

thermosyphoning is predicted for the less probable accident 

scenarios only. 

Thermosyphoning requires a heat sink at 
the boilers. Water is provided either by the normal 
feedwater system, the auxiliary feedwater pumps or, on 
de pressur isat ion of the secondary side, by an emergency 

supply system. 

1 .5 Review of past work 

The problem of heat transfer and 
fluid flow in natural convection loops h; stimulated 

considerable interest in the past due to its important 
applications in oceanographic and geophysical phenomena as 
well as in such practical systems as nuclear reactors and 
solar heaters. The early studies on such systems started 

with Rayleigh's analysis concerned with the stability of a 
fluid at rest and the onset of motion due to a 
temperature gradient in the direction of the body force. A 
review of external flows caused by buoyancy forces was 
presented in C23. Another type of free convection flows is 
encountered in enclosures and cavities the review of which 
is given in C33 and C4U.An overall study of the 


buoyancy- 

-induced 

flow and 

transport 

is 

given 

in the text 

C5U for 

various 

physical 

cases. 







Some of 

the 

previous studies of 

natural 

circulation loops 

were concerned 

with 

the stability 

of the 

steady- 

-state motion .Kel ler 

C63 

and 

Welander C7D 


considered a simple rectangular thermosyphon consisting of 



two insulated vertical branches with point heat source and 

sink at the centre of the lower and upper horizontal 
segments respectively.lt was pointed out in these 
pioneering works that, when a certain parameter representing 
the ratio of the buoyancy force to the frictional force 
exceeds a critical value, an oscillatory motion of the 
fluid may occur. The work of Creveling et al . considered a 
toroidal loop and demonstrated experimentally and 

theoretically the presence of instabi 1 it ies .Studies of 
stability behaviour of circular loops have also been done 

by Zvirin et al C9U.Zvirin et al. C1Q3 have also studied 
the steady state, transient behaviour of a natural circulation 
loop with two vertical branches with point heat source 
and sink and have shown that the stability depends on 

the type of the temperature distribution assumed. 

The loop configuration is one 
of the few alterable factors in thermosyphon design. ft 

detailed study of this factor is done by ChenCIIU for 

better understanding of the operational and safety problems 
caused by the periodic instability in the natural 
convection loops. 

The instability associated with the 

onset of motion in a thermosyphon has been studied by 

ZvirinC123 to show that the rest state is always unstable 

as the Rayleigh number exceeds some critical value. 

The above studies considered one 

dimensional analysis of the thermosyphon with the 
coordinate running along the loop.liertol et al.C13H tried a 

two-dimensional analysis of the toroidal thermosyphon and 

Mallinson et al.C141 analysed the three-dimensional case of 
a rectangular thermosyphon. 

HartC153 gave a new approach to the 

thermosyphon analysis by solving the conservation equations 
by a set of coupled master and slave ordinary differential 

equations and by making the chaotic analysis, he claimed 

the instabilities to be chaotic rather than oscillatory. 

The effect of viscous dissipation 

is studied by Zvi rinC16U .The transient, steady state and 
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stability behaviour of a toroidal the rmosyphon with 
throughflow is studied by Mertol et al.C173. 

As far as the application of 

thermosyphon to decay heat removal in nuclear reactors is 
concerned, the text by LewisC18U threw some light for the 

different formulations .The text by Guppy et alC193 discussed the 
decay heat removal in fast breeder reactors by natural 

convection in detai 1 . Zvi r inC203 reviewed the transient , steady 
state and stability behaviour of different natural 
circulation loops including PWRs. 

Ghosh et alC21Z) studied the 

natural circulation cooling of PHWR for decay heat removal 
after a reactor trip under loss of class IV power. They 
showed that under bottled condition, thermosyphon cooling is 
adequate to remove upto 10% of full power without boiling 
and upto 12% with boiling in the primary coolant 
channels . 

Spinks et al.C22T considered the 

thermosyphoning in CANDU reactor .ChatoC233 studied the 

natural convection flows in parallel channel systems . Zvi r in 

et al.C24T investigated a natural circulation sytem with 
parallel loops both analytically and exper imental ly .They 
studied the steady state and transient characteristics and 
the effects of core flow resistance, power distribution and 
upper plenum design. Flow oscillations were observed under 

certain conditions which were accompanied by instabilities 

and flow reversals. 

The transition from forced flow to 
natural circulation was studied experimentally by Gillette 

et alC25U.They performed their experiment to simulate 

EBR-II,a liquid-metal-cooled fast breeder reactor. The test 

was initiated by abruptly tripping an electromagnetic pump 

which supplied 5-6% of the normal full operational primary 

flow rate. The flow coast-down reached a minimum after 

which the flow increased as natural circulation was 
establ ished.The effect of secondary system flow through 

the IHX and reactor decay power level on minimum in-core; 
flow rates and maximum in-core temperature were examined. 


7 



1 .6 Scope of present work 

In the present work, a nuclear 

reactor primary heat transport system has been simulated 
as a rectangular loop in a vertical plane the bottom 
segment of which is heated and the top horizontal segment 
is cooled by passing a coolant. In the present analysis, the 

approach is made initially by taking simplified assumptions 
and latter on the complexities being handled by relaxing 

the simple assumptions. 

In Chapter-2, to get the feel of the 
comlex nature of the problem as regards to stability, the 

simple case of a rectangular thermosyphon with point heat 
source and point heat sink has been analysed in 
detail. Two approaches have been made. In the first case, a 
linear temperature distribution is assumed which fails to 
give any information regarding the stability though steady 
state and transient behaviours are predicted correctly.In 
the second approach,a numerical solution of the equations 

show the in s tab i 1 ity associated with the problem. A linear 

stability analysis is also done. 

Chapter-3 deals with the real 

situation of extended heat source and sink. For 

simplicity, the sink is assumed to be at constant 
temperature whereas three different types of heat source 
are chosen to give uniform, triangular and sinosidual heat 
fluxes. The analysis is done considering one-dimensional flow 
with the coordinate along the loop. The fluid flow and 

heat transfer character istics have been analysed for each 

case . 


more practical 
presence of a 
of a constant 
previous chapte 


Chapter-4 takes 
situation by analysing 
heat exchanger as the 
temperature heat sink 
rs . 


into account the 


the 

problem 

with the 

heat 

sink 

in stead 

analysed 

in the 



Experimentally the flow has been observed to be 
two-dimensional .Again for one-dimensional analysis, we have to use 
some correlations for the friction factor and heat transfer 
coefficient. As the appropriate correlations for the natural 
circulation flows hare not been developed, we have to use the 
correlations available in forced convection flows for those 
quanta, t* es. Thus the combined effect of these two defects in our 
ea; lies* ccd<_ls may be the cause for the discrepancies between the 
Lnu-ot at* c-1 and the experimental results. So in Chapter -5, a two 
dimensional analysis of the problem is done. 

F'.ncll> Chapter -6 deuls with the future scope of the work and 
draws some conclusions of the present work. 
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CHAPTER-2 

RECTANGULAR THERMOSYPHON WITH 
POINT SOURCE AND SINK 


2.1 Introduction 

Before analysing the real problem, it is 
advantageous to consider a simple model to have a better 
insight into the mechanism of convection. The model is 
materialised by taking a narrow tube of uniform 
cross-section and forming it into a closed loop. The tube 
is filled with fluid that is kept well mixed over the 
cross-section .The fluid is subjected to given heat sources 
and sinks along its path. The fluid motion is driven by 


buoyancy force 

and 

will also 

be retarded by 

frictional 

forces .The 

tube 

is 

considered 

to be insulated 

everywhere 

except at 

the 

top 

and bottom 

,where over short 

distances s 

temperature 

of 

the 

wall is 

kept at -AT 

and AT 

respectively 

-Thus 

we 

conside r 

the limiting 

case where 

As >0 while 

k 

>o> in such 

a way that the 

heat flux 


remains f inite, i .e, point source and point sink. 

In this chapter, the steady 

state, transient and stability behaviour of a system shown 

in Figure-2.1 as comprising of two vert: si branches with 

point heat source and sink has been discussed. 

2.2 Formulation 

Let us consider a portion of the 
fluid with uniform cross-sectional area A and length L as 
shown in Figure-2. 2. The fluid is driven by the pressure 

difference between the end points and by a buoyancy force 
and is retarded by frictional force. The following 
assumptions are made. 

DThe fluid is incompressible due to the small pressure 

change encountered in natural convection. 

2) The flow is one-dimensional and the coordinate runs 

along the loop. 

3) A11 the fluid properties and heat transfer coefficients 
are constant. 

4>Boussinesq approximat ion is val id, i . e ., density is assumed 
to be constant in all the conservative equations except 


10 



-^7 




^IGuftE Z-l-The Simplified 

model 


'X 

t 



FlGUftE • Tube segment 



for the body force term, caused due to density 
difference, where linear variation of density is assumed. 


i . e ,p- p ZA-ft ( T— T >3 (2.2.1) 

o o 


5) The 

fluid is 

well 

mixed 

over 

the 

cross-section. The 

characteristic 

time 

for 

mixing 

is small 

compared to the 

time 

required 

for 

the 

fluid 

to be 

advected an 


appreciable distance along the loop. 

6>Due to well-mixing ,the temperature is assumed uniform 
over the cross-section. 

7)The tangential friction force is a function of the 

instantaneous flow rate. Here laminar flow is studied. Hence , a 
linear relation is assumed. 

8>The heat flux between the tube and the fluid is 

proportional to the temperature difference between wall 
temperature and fluid temperature T. 

Considering the tube segme n t in 

Figure-2.2, the following conservation equations are derived. 


Continuity: V.v =0 (2.2.2) 

Momentum: A 

——p g e -V p+ F (2.2.3) 

Energy: 

== =k CT ( s )— T3 (2.2.4) 

dt o 


This constitute a set of 3 
equations with 4 unknowns v,T ,p and p.So this equation 
set along with eq. (2.2.1) relating density and temperature 
constitute the total set of equations. 

For one-dimensional case, eq. (2.2.2) 


gives : s v 

___ = o ,i.e., v = v(t) (2.2.5) 

S $, 

So velocity is uniform along the loop. 
Taking component of (2.2.3) along the loop. 


P <4 V ~~pg cos 4> - S p /d s + F (2.2.6) 

at a 

Integrating over the fluid volume. 


0 d 

o ^ SSv dA ds =d\f-p gcos$> dA ds — ff (d p/d s ) dA ds+// F^dAds 

As flow rate Q=Sv dA is constant by (2.2.5)? ds cos^-dz 
and using (2.2.1) and assumptions (6) & (7), we get. 


il 



(2.2.7) 


p o L dz+cr dz-AT dp-p LRQ 

Ut O w 

where R is a frictional coefficient 

Now for a closed tube, eq. (2.2., 7> reduces to 


p L(dG/dt)=p A9g#T dz- p LRQ 
o o o 

as # dp=0 and # dz=0 

i.e., dQ/dt=A/?g/L # T dz - RQ (2.2.8) 

The energy equation gives 

id 1/d t)+Q/A id 1/d s)=k(T -T) (2.2.9) 

O 


The tube is symmetrical with respect to 
the vertical. The temperature is antisymmetric with respect to the 
center of the loop, i . e . ,T( s+L/2)=-T(s ) .These conditions are shown 
in Figure-2.3. 

Applying the symmetry and antisymmetry 
conditions to the integral tf T dz,we get, 

. _ . _ ,L /2 • 

# T dz =2 S T ds 

Thus equation of motion becomes 

dG/dl=2Ag/5/L /^“t ds-RQ (2.2.10) 

The solution of (2.2.9) and (2.2.10) need some boundary conditions 
involving temp. The temperature of the fluid coming out of the 
source or sink depend on the temperature of the ingoing fluid and 
flow rate .Appendix-! proves that if a fluid particle passes through 
the source or sink for a time At, then the difference in 
temperature between the outgoing and incoming fluid is given by 
T -T. = (AT— T. ) (1-exp (At)) 

out in «.r» 


i.e.,T -T = (AT— T > ( 1-exp(-kA sA/Q) ) for source (2.2.11) 

out in ih 

= (— AT-T. ) ( 1-exp( kAsA/G ) ) for sink 

vrt 

For the computation in the range 0<s<L/2,we need to know T at 

in 

source when Q>0 and at the sink when Q<0. 

By antisymmetry ,T. =T(L-As/2)=-T( (L-As > /2) f or first case 

ir» 

=T( (L+As)/2)=— T(As/2) for second case 
In the limit as As >0,T. =— T(L/2) for first case 

ITS 

=-T(0) for second case 


So , eq . (2.2. 1 1 ) reduces to 

T +T =(AT+T ) ( 1-exp(-kA!kS/G> ) -,G>0 

= (AT+T ) ( 1— exp(-kAAs/Q) ) ,G<0 
•so r 


j (2.2.12) 
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T(s+iQ 

= -TC*) 





S = 4^ s 


Figure 2-3 Coordinate alon^ 
tbe tube and antisymmetry 
Condition 


£ 


♦5 


4 
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Fjgube a -4 Discretisation 
of the loo|? 



£.3 Non-dimensional i sat ion 

The equations are non-dimens ional i sed by 

the transf or mat ion 

s*=L/2*s jt* = (L/(2kAs) >t ;q*=k£Asq ;T 4T#T 

and finally take the form 
dQ/dt+«Q=«x^ A T ds (£.3.1) 

O 

a t /a t +q a T/a s =o <£.3.£> 

T +T = ( 1+T ) (1— exp(— 1/Q)) . QX) (£.3.3) 

•«0 ' I 

=(— 1+T ) ( 1-exp( — 1/G) ) , GU .0 

•=o 

where the non-dimensional parameters are given by 

2 

«=Gravity parameter= g^ATL/ ( kAs ) 

^=Friction parameter=RL/ (2kAs ) J (t.3.4) 

£.4 Steady state solution 

We assume that Q > O, hence the motion is 
upward in O < s < I.For steady state solution, the time derivative 
vanishes . 


Eq.(£.3.£) reduces to T=constant=T (£.4.1) 

Eq .(£.3.1) results in ^Q=o*T (2.4.2) 

Eq. (£.3.3) gives 2T=(1+T> ( 1-exp(-1/G) > (£.4.3) 

Eliminating T between (£.4.2) and (2.4.3) , we , get , 
£G/(G+a/£r) = 1-exp(-1/G> (2.4.4) 


The LHS and RHS of (£.4.4) show that LHS is a monotoni cal ly 
decreasing and RHS is a monotonically increasing function of 1/G. 
LHS decreases from 2 to 0 and RHS increases from O to 1 as G 
decreases from oo to O.So there is a point where they intersect 
giving rise to steady state value G. 

Eq. (2.4.4) also show that the steady state 
solution is a function of a/e . Table-1 gives different values of Q 
for different a/e ratios. 

£.5 T ransient solution 

Eq. (2.3.1) and (2.3.2) form a set of 
coupled integro-dif f erential equations in G and T.The temperature 
being a function of both position and time needs a boundary 
condition in space given by Eq. (2.3.3) .So to have a transient 
solution, we need initial conditions only. Here we consider the 
case where the motion starts from rest. 

So at t=0,G=0 and T= 1 at s=0 

=0 at all other s 
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a/e 

9 

T 

0-100 

0-100 

1-000 

0-500 

0-417 

0-834 

1-000 

0-648 

0-64S 

2-000 

0-95S 

0-479 

10-000 

2-21S 

0*222 


Table 1 . Values of non-dimensional flow rate and 
temperature in steady motion 


The transient behaviour has been studied by two 

methods . 

DZvirin’s ordinary differential equation method 
2>Numerical solution method by finite difference 
2.5.1 Zvirin's method 

This method comprises two stages .During the 
first stage, a linear temperature distribution is taken from the 
bottom to a penetration height H(t) and T=0 for H < s < I.The 
penetration height increases with time and the first stage ends 

4 it 

when H(t)=1 at some time t=t .Thus the second stage begins for 
t > t when a linear temperature distribution is assumed over the 
whole range 0 < s < 1. 

Reduction of governing eg . into ODEs s 

Integration of energy eq. (2.3.2) gives 


i / 

dt o 


T ds -GT =0 
o 


Assuming linear temperature var iat ion , i . e . ,T=T — ' T s/H ,we get 


o 4* < HT > = 

2 dt o 


Momentum equation 


QT (2.5.1) 

o 

(2.3.1) reduces to 


— + cG =1/2*otT H (2.5.2) 

dt o 

The boundary condition (2.3.3) becomes 

T = 1-exp(-1/G) (2.5.3) 

o 

The above equation set represents a set of ordinary 
coupled differential equations with 3 variables Q,T^ and H with 
initial conditions G=0 , T =1 , H=H at t=0 (2.5.4) 

O i 

where H is a small value representing the initial asymmetry of 
heat source. 

The above set can be reduced to a single equation by 

dif f erientiating(2.5.3)with respect to t and substituting 

expression for T from (2.5.3) .This results in 

o 

“f tb +s =<xQ( 1-exp(-1/Q) ) (2.5.5) 

d t dt 

with initial conditions G=0 ,dG/dt=1/2 otHat t=0 — (2.5.6) 

For the second stage, i.e.,t > t* 
integration of (2.3.2) gives 
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■7 S* T ds -Q(T -T ) =0 <2.5.7) 

dt o to 

With linear temperature distribution,./* 4 T ds=1/2(T +T >=T 

r O O i m 

Now using <2.3.3),T =1-2T exp(-1/Q> / < 1-exp<-1/G) > 3 

O w / o 

T =2T /<1-exp(-1/Q) ) — 1 J 

4 m 


Now (2.5.7) reduces to 

dT /dt=2G C1-(1+exp(-1/Q) )*T / ( 1-exp<-1/G) ) 1 <2.5.9) 

m m 

The momentum eq<2.3.1> reduces to 
dQ/dt+£:Q=otT <2.5.10) 

TO 

These two coupled ordinary differential eq. in G and T are to be 

TO 

solved with the initial condition 
Q=G ,T =1/2T at t=t* 

m L cl Of* — X to O 

Solution Technique? 

The non-linearity of equation (2.5.5) in stage-1 and 

eq.s <2.5.9) and (2.5.10) in stage-2 eliminates any chance of 

analytical solution. So they are to be solved numerical ly. As the 

equations are ordinary differential equations and are of the 

initial value type, they are solved by Runge-Kutta method. 

The solution algorithm first solves 

the equation for the 1st stage -Start ing from t=0,at each time 

step,Q,dG/dt are calculated. From<2. 5. 3) T is found out and 

o 

from<2.5.2) H is found. The first stage solution continues till 

* * 

H=1,the corresponding time being t .Taking the value of Q at t=t 
as initial value, the coupled set(2.5.9> and (2.5.10) are solved at 
each time for G and T .From (2.5.8), T and T are found. The process 

m 0 4 

continues till steady state is achieved. 

Stabi 1 i tv Characteristics 

Let us consider small deviations from steady 


state in the form 

Q ( t ) =G+Q * ( t ) T (t)=T +T ’ <t) (2.5.11) 

TO TO TO 


where T and G are steady state values given in section 2.3. 


TO 

Introducing (2.5.11) into (2.5.9>and 

<2.5. 10) , subtracting the steady state relations and using the 
linearized stability procedure, we get, 

< 1-m) dT Vdt= -2G CmG’/Q 2 ( 1+T > + ( 1+m)T *3 <2.5.12) 

m mm 

dG*/dt=— eG’+otT ' (2.5.13) 

m 

where m=exp(-1/G> 



To reduce the above equations into algebraic ones and to check 
their stability, we assume. 


T ’ =T exp(o-t) G ’ =Q exp(O't) (2.5.14) 

mm 

where o', in general, is a complex number. 

Now, we get,(1-m)o-T = -2Q CmQ/Q 2 ( 1+T >+(1+m)T 3 (2.5.15) 

fft Tfi Th 

&Q=-£Q+ aT (2.5.16) 

m a 


Elimination of T and Q gives rise to the 

m 

characteristic eq. 

o-*+o-(*+2G/T )+2cCG/T +2m/ ( 1-m)* 3=0 (2.5.17) 

m m 

As the coefficients in the quadratic 
equation are always pos i t i ve , the re are no roots of o'with positive 
real part. Hence the solution is stable. 

Results and Discussions 

The plot of heated height H,flow rate Q,T 0 and 
as a function of time has been shown for various values of the 
parameters for different initial penetration depths in Graph 

1.1 to 1.8. A close analysis of the results show the followings 

(1) The time required for the full growth of the heated region of 
the loop is more for low initial penetration depth as expected. For 
high initial penetration depth , i . e ., wi th some initial flow, the 
natural circulation flow establishes faster. 

(2) As the ratio of gravity to friction parameter is increased, f low 
is established at much less time. This is due to the increase in 
the driving buoyancy force and decrease in the retarding friction 
force. 

(3) The study with linear temperature distribution assumption : 
explains steady state and transient behaviour well but fails to | 
explain the instability associated with the flow. Thus it may be 
concluded that the stability of such a system depends on the j 
assumed temperature distribution. 

2.5.2;Numeri cal solution by Finite Difference J_ | 

Zvirin’s method failed from certain drawbacks. j 
(1)The assumption of linear temperature profile is much a 
simplified assumption .The temperature distribution may be linear ! 
at the initial period when convection is not significant and the 
mode of heat transfer is conduction. But as natural circulation 
flow becomes significant, temperature distribution will not be ; 
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simply linear. 

(2)The linear temperature profile assumption fails to explain the 
instability associated with the process as shown earlier. 

In order to get a true picture as regards 
to stability, it is better to take up the original equations and 
make a numerical solution of these equations though it may require 
some additional mathematical analysis. 

Solution Technique 

The coupled non-linear integro-differential 
equation set(2.3.1) and (2.3.2) involving G and T are solved by 
finite-difference scheme with the boundary condition and initial 
conditions as prescribed earlier. The discretisation of the loop is 
shown in Figure-2.4. 

The integral in eq. (2.3.1) is 
decomposed by Simpson’s rule. An implicit scheme is used for the 
time deri vat i ve .The upwind scheme is used for the space 
derivative, i .e, backward difference scheme in the present case. 

Thus the finite— difference form of the eq. (2.3. 1 )and 
(2.3.2) are 

N 

, _ t •+■ 4 _ t . ... _ t ■+ 4 . , _ _ _ t+4 

(G -Q )/At+.eG =c*As/3 T a. T. 

V. t 

where a. =1 for i=1 and N 

i. 

=4 for even i and 2 for odd i 


. _t 4- 4 t _ t + 4 ,_t+4 _t-*-4 . , . _ _ _ _ . , 

(T. — T. )/A*t+Q (T. —TT )/As =0 ; for i=2 ,3 r » • . »N 

L X. V L— 4 

A _t -HI . . _ .t+A . , . , . ,, i .■+•! . . 

T +T =( 1+T )(1-exp(-1/Q )) 

AN N r 

N 

i.e. ,G t+1 = (G t +c*AsAt/3 T a.T t+1 )/(1+cAt) (2.5.18) 

i=1 1 1 


t + i t t + 1 t + i t + i 

T =(T +AtG T /As)/(1+G At/As ) ; i=2,3, . .N (2.5.19) 

\. t v - 4 

T^ +1 =1-( 1+T^ +1 )exp(-1/G t+1 ) (2.5.20) 

1 1 

For solution, first eq. set (2.5.19) is 

i X 

solved. The non-linearity arising is handled by taking Q =Q and 
T^ =T^ .After finding T^ ,....T^ , corrected value of T^ 

is found from(2.5.20) and then corrected value of G t+i is found from 
(2. 5. 18). The process is repeated till convergence .The marching in 
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time is done till steady state is achieved or till a particular 
time to see the flow oscillations if there is no steady state. 
Results and Discussions 

The results are given as the plot of flow rate 
vs. time, temperature at various nodes vs. time and temperature 
distribution along the loop at different times for several sets of 
parameters a and e in Graph 1.9 to T. 20. The flow rate curve shows 
monotonous increase till steady state for <*=0.4 and «=0.2 whereas 
for other ranges <x=2,e = 1 j<x=20,e=3 and <*=40,o=6 flow rate comes to 
steady state after some initial osci 1 lat ions .So for some range of 
<* and <,we may encounter flow oscillations or even flow reversal. 

The temperature distribution vs. 

time shows that near the source, temperature decreases gradually to 
steady state value whereas at other points temperature increases 
till a maximum and then falls to steady state value. This can be 
explained as follows. At initial time, when flow rate is 
smal 1 , i . e . , natural convection has just been establ ished , the 
temperature of the outgoing fluid from the source becomes close to 
the source temperature as the fluid takes more time to pass 
through it. As flow rate increases at subsequent time, the outgoing 
fluid temperature is nearer to the incoming fluid temperature which 
is relatively cold. 

For other points , in it ial ly conduction is 
predominant mode of heat transfer. As convection is establ ished, the 
temperature rises due to the advected hot fluid particles emanating 
from source till a time after which the temperature decreases due 
to relatively cold particles coming out of the source , mixing with 
the hot particles and being conve cted. Final ly steady state is 
achieved when temperature throughout the loop is constant. 

2.6 Stability Characteristics 

The fluid motion is driven by buoyancy force and 
retarded by friction force. The buoyancy force integrated around the 
loop is B=p Q gA# T dz (2.6.1) 

The temperature distribution of the fluid is solved from the 

energy equation Q/A O T/6 s=k CT^(s)-TII (2.6.2) 

where T^ ( s ) is prescribed wall temperature. 

So temperature distribution is a function of G and 
hence B also. To find the behaviour of B with G,we f ind,f rom(2.6.2) 
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as Q — >w t S7/Ss~ — >0,i.e.,T is uniform and hence B — >0 

For G=0,T=T and hence B is non-zero if # T dz is not zeroCCase-1) 
o o 

and B=0 if# T dz=0(Case-2) 
o 

The frictional force F is a function of G 
and increases with Q for any normal friction law. This is also 
shown in Figure 2.5. 

Now, clearly in both cases, the two curves 
must meet at least in one point. In case-2, zero motion is also a 
possible steady solution. 

In Figure-2. 5a, as G is increased above its 
balanced value, buoyancy force decreases and friction increases 
thereby counteracting the changes. Thus the system seems to be 
stable . 

In figure-2. 5b, if rate of increase of B is larger than 
that of F near G=0,i.e.,if dF/dQ < dB/dG at G=0,then increase of G 
results in larger increase of B than that of F. Hence instability 
may result. But as Q is increased further, we get a point after 
which steady state prevai Is . Thus it seems that the system would 
always end up in a stable situation. But experimentally instability 
has been observed when buoancy forces are large. This is explained 
analytically in the following section by linear stability 
analysis . 

Linear Stabi 1 itv Analysis 

Let us analyse the equation set 

(2.3. 1 )-(2.3.3) .By taking a first order perturbation, 

Q(t>=Q-K5G* (t) T ( t , s ) =T ( s ) +6 T * (t,s) (2.6.3) 

and applying to (2.3. 1 ) — (2.3.3) , we get, by sepatrofclon of zero and 
first order terms, 

*G=ot T ds (2. 6. 4. a) 

o 


Q ST/ds=0 (2.6.4. b) 

(2.6.4. c) 

(2. 6. 5. a) 

dT* /dt+Q'0T'/ds=O (2.6.5. b) 


T +T =(1+T ) ( 1-exp(-1/G) )• 

fi=u e= i. 


dGVdt+aG*=oi f T' ds- 
o 
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T * +mT ’ +nG'=0 (2.6.5. c) 

•sO *5=4 

where m=exp(-1/Q ) = ( 1-T) / ( 1+T ) ; 

n=(1+T)exp(-1/Q)/(Q*Q)=(l-f )/(G*G> (2.6.5. d) 

Equation (2.6.4) represents the steady state case whose solution 
is given in section (2. 4). To solve (2. 6. 5), we take, 

A A 

Q " ( t ) exp (<>t ) T * ( s , t )=T ( s ) exp(o , t ) (2.6.6) 

where o', the stability parameter, in general, is a complex number. 
Applying in (2. 6. 5), we get. 


(<y+t ) Q=ec/‘ i T ds (2. 6. 7. a) 

o 

<yT+G dT/ds=0 (2.6.7. b) 

T(0)+mT(1 )+nQ=0 (2.6.7. c) 


By solving for T and Q, 

T=C exp(-<>s/Q) 

Q=o*/(o‘+£t) ds^CQa/Co 1 (<y+E ) II ( 1-exp(-tJ'/Q) ) 

o 

Now B.C. (2.6.7. c) becomes the characteristic equation 

1+m exp(-' 0 '/G)+naG/ (o- {&+£ ) ) E 1-ex p( -o-VG) 3=0 (2.6.8) 

For real , posit ive o^all terms in LHS of (2.6.8) become 
positive and (2.6.8) does not hold. Hence there is no possibility 
of real , positive o 1 , hence no instability of exponential growth 
type. But there may be oscillatory instabi 1 i ty .The neutral 

stability can be analysed by taking as a purely imaginary 

quantity iw. 

j\i *\j __ 

Taking «=not/G ;s=s/G ;u=w/Q,we get, 

<V 'V < \r ■'Vr 

( e x p ( i<*> ) +rri )/(exp(i^) — 1> / ( io* ( > > =0 

*v 

Separation of real and imaginary parts alongwith expressing otin 

*Yr 

terms of T and e only and substituting for m from (2. 6. 5. d), we 


get, w 2 + (c— A/2) Z = (A/2) 2 (2.6.9) 

s=-1/T Z cot (Z/ 2) (2.6.10) 

* 

where A=A(T)=-( 1-T* )/T 2 InC ( 1-T ) / ( 1+T) 2 (2.6.11) 


The plot of A vs. T and —u> cot(w/2)is shown 
in Figure-2.6 and 2. 7. The curve represented by (2.6.9) is a circle 
and by (2.6. 10) is a multi-branched curve, both being functions of 
T in Z -s plane. Thus for certain T=0. 265, both curves first touch 
each other.lt is the largest T for which flow becomes unstable 





corresponding to smallest possible flow rate G=1.84 and angular 
frequency of os cillation 3.6G.As angular frequency of fluid 
particle is riQ, thermal disturbance travels slightly faster than 
the fluid. 

As Q is increased, there are two intersection 
points of the curves of (2.6.9) and (2.6. 10) .Thus there are two 
values of £ for which neutral oscillations occur if the value of Q 
is above the critical value;in between these ranges the solutions 
are amplif ied.The neutral stability curve is shown in Figure-2. S 
in plane. 

Physical explanation of instabi 1 itv 

The instability can be explained physically 
by considering the interplay of buoyancy and friction forces. The 
instability arises as these two forces may not be in phase at each 
time . 

Let us consider the system to be operating in 
steady state. If a small thermal disturbance causes a pocket of 
fluid to emerge from the heated section slightly hotter than is 
normal for this steady state, the buoyany is momentarily increased 
and the flow accelerates slightly. This hotter-than-normal 
pocket proceeds through the cooled section more rapidly than normal 
and is warmer than normal as it leaves the sink. The total buoyancy 
force is therefore decreased decelerating the flow. Thus the pocket 
enters the heated section once more, now at a higher inlet 
temperature and a lower velocity. These two effects combined 
together cause the pocket to emerge from the source at a higher 
temperature than before. The same logic applies for a pocket of 
fluid which initially emerges from the source cooler than normal 
leading to the same oscillatory result. 
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CHAPTER -3 

RECTANGULAR THERMQSYPHQN WITH 
EXTENDED SOURCE AND SINK 


3. ^Introduction 

The analysis of the previous chapter gives some 
insight into the complexity of the problem. However , the idea of a 
point source and point sink is not a real one, only assumed for 
mathematical simplicity in analysis. The real problem involves a 
line source and a line sink. So in this chapter, we consider a 
rectangular thermosyphon loop in vertical plane with insulated 
vertical branches whose lower segment is being heated and upper 
segment cooled. Here three types of heat fluxes are 
considered, i .e. , uniform, triangular, sinosidual .The heat sink has 
been assumed to be at constant temperature for simplicity. 
3.2:Formulation 

The loop to be analysed is shown in Figure-3. 1 .The 
following assumptions are made. 

(1) Due to small pressure changes in natural convection, the fluid 
is incompressible. 

(2) The flow is one-dimensional with the coordinate running around 
the loop. 

(3) The fluid properties and heat transfer coefficients are 
constant . 

(4) Boussinesq approximation is val id, i . e , densi ty is assumed to be 

constant in all conservation equations except for the body force 
term in the momentum equation where linear variation of density 
with temperature is assumed, i .e ,p= p C1-Y?(T— T >3 (3.2.1) 

(.?. 'he fluid is well mixed over the cross-section. 

(6)The temperature is uniform over the cross-section due to 
mixing . 

(7>The frictional force is a function of instantaneous flow rate. 

Now, the governing equations are: 

Continuity: d v/tfs=0, i.e,v=v(t) (3.2.2) 

Velocity is uniform around the loop. 

Momentum: dv/dt= - pg cos3-dp/ds-4T^ /d 

where <9=Angle between vertical and tangent to the loop at a point. 
Now to eliminate tfce pressure term, integration around the loop is 
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done. Thus, we get, 

p dG/dt # ds/A(s)=— g# pdz — t 4t /d ds 

V V 

Now using the relation (3.2.1) and using t =f/2 p G 2 /A 2 , 

V V 

we get, p^dQ/dt #ds/A(s)=gp^/9# T dz— #2fp Q 2 /(dA 2 ) ds 
Now assuming uniform cross-sectional area around the loop 
p dG/dt L/A=p gfi* T dz-p /2RQ 2 

where R=4tf f / ( dA a ) ds=Overal 1 friction factor, which accounts for 
both friction and form losses 

Thus ,L/A dG/dt=g/j# T dz-RQ Z /2 (3.2.3) 

Energy: Neglecting viscous dissipation and axial conduction, 

P c ~)=-3 r o <s <W/2(Unif orm heat flux) 

v p H v 5 a 


=0 , W/2 <s< ( W+H) /2 and W+H/2<s<W+H (3.2.4) 

4h (T-Tv ) 

L “ , ( W+H ) /2< s< W+H/2 
d 

Eq. (3.2.3) and (3.2.4) represent a set of coupled 
integro-dif f erential equations in two variables G,a function of t 
only and T,a function of s and t. Hence we need boundary condition 
and initial condition . 

B.C. Continuity of T around the loop 


(3.2.5) 


I.C.sAt time t=0,Q=0 and T=T for all s 

o 

3.3 :5teadv state solution 

The steady state equations are derived by making 

time derivative zero. So we get f rom(3.2.3) , (3.2.4)and(3.2.5) 

Momentum: g/9* T dz=RQ 2 /2 (3.3.1) 

Energy: p c ? |-=-§ ,0<s<W/2(Sour ce ) 

w p A as a 

=0 , W/2< s< ( W+H ) /2 and W+H/2<s<W+H( Vert i cal side) 

, (W+H) /2<s<W+H/2( Sink) (3.3.2) 

d 

B.C . :Continuity of T along the loop 

To solve the above set of equations, we first 

solve the energy equation to find temperature distribution with Q 

as a parameter and then apply this distribution in (3.3.1) to find 

steady state values. 

Solution of energy eq uation: 

Source (0<s<W/2) :T=--^~~-+C 

p C dG 1 

v p 

U LJ+H H 

Vertical branches (=<s<-=-,W+=;<s<W+H) :T-constant , say 

4 c lj U-fH 

=T H for g<s<-g-(Hot segment) 

=T for W+g<s<W+H(Cold segment) 

Sink(«|B<s<W + |) =T=C 2 exp<-^§S 55 )+T v 


. It 



To find the constants,we apply the boundary 


condition. 

Continuity at s=W/2 : “§~dG2 +C 1 =T h (3.3.3) 

v p 

a t C -W±H , r o . 4hA(W+H) . _ T 

at s- 2 .C 2 exp( ------- ) +T v -T h (3.3.4) 

w p 


, H 


4hA(W+ l ^> 


at s=W+- s T c =C 2 exp(~~-~2--)+T , (3.3.5) 


P C Gd 

w p 


V 


at s=W+H: T =C - 

c 1 


(3.3.6) 


Eq. (3.3.3) and (3.3.6) give T -T 
^ a h c 2p C dG 2p C G 

V p y p 


Since , power=P=qndW/2 r AT =Temperature difference across heated 

R P 

sect ion=T -T = — --- <3.3.6) 

h c p C Q 

w p 

Eq. (3.3.4) gives C =(T -T )exp(~-~”-) 

M 3 2hv k 2pCG 

w p 

Eq. (3.3.5) gives T =(T -T ) ex p 5 +T (3.3.7) 

3 c h v 2p C G v 

w p 


Now temperature distribution is given by 

T=— |-“+T (Source) 

P C WQ c 

w p 


=T (Hot 

H 

= (T -T 


H V 


segment ) 


) exp( 


hndi«+H2 

Zp C G 

w p 


)+T 


(Sink) 


(3.3.8) 


=T c (Cold segment) 

Now momentum eq. (3.3.1) gives 

g/9 (T -T )H/2=RG Z /2 (3.3.9) 

Eq . (3.3.6) , (3.3.7) and (3.3.9) are to be solved for unknowns T r T 

H C 

and Q. 

From(3.3.6) , f--“ ? r= R G 2/ 2 

C 8 

v p 

So,Q=(^““-)^ 3 (3.3.10) 

Vr K 
v p 

Now T and T can be solved and temperature distribution in the 

H C 

loop can be found from (3.3.8). 
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3.4 Non-dimensionalisations 


* I=I W 

1 “q/h~ 


The variables are non-dimensionalised as following 


Now 


* _ s 
s H/2 
2 

dG_ Qch_ dGL 

dt ( H+W ) A dt 


q*=-Q— 

Qch 


t *__tQch_ 
(H+W) A 


£I = Qch3 ST* 
St hLA St 


*T_2g ST* 
Ss hH as 


ds=H/2 ds 


Substituting these relations in (3.2.3)— (3.2.5) , by dropping * from 

variables and defining aspect ratio =B=W/H,we get, 

. dQ __ i + B_ , 2 + 2B , _ „2 

1/F j 7 = GC/ T ds -S T ds3-G 

dt b *+a b 


2‘(1^Bl fl +Q |J=' y ',0<s<B(Source) 

=0,B<s<1+B(Hot segment) 


=— MT , 1+B<s< 1+2B ( Sink ) 

=0, 1+2B< s<2+2B ( Co 1 d segment) 
with continuity of T along the loop and 
initial condition of Q=0 and T=0 at t=Q 
with the non-dimensional parameters defined by 
F=Friction parameter=RA 2 /2 

G=Modif ied Grashof numbe 

RhG ,ndB RhG, 

Ori Ori 


M=Geometry parameter= 


2hHA_ 
PC Q d 

p ch 


(3.4.1) 


(3.4.2) 


In natural convection , the main problem encountered in 
non-dimensional isation is the non-availability of a character isti c 
flow rate as compared to forced convection .Hence this problem is 
solved by taking the value of G to be equal to the steady state 

c h 

flow rate. Hence, G , =Q 

ch am 

To find the non-dimensional parameters,we 
need to know R,G and h.From definitions 

c h 

R =-™2 where f=Friction factor, a function of Reynold’s No. 
dA 


a , 

b b 

A 

Re 

(Q u pdV 

ch 

' * R : 

4AUi b 

7 1 v 

2 b b. A + b 

A Q (£> d 

ch 


(3.4.3) 


The relations (3.3.10) and(3.4.3> involve two unknowns R and G cK . 


2/5 



2-b i+b 

By solving, Q c h = ^ 3^—-------- — ) 4 ' <a-b> (3.4.4) 

4C Laju p 

p 


For laminar flow r a=16, b=1, 

Q = { s£HP0dl YJ(1 21"™- > 

'■* I ' / * « i t “ \ r-ir— / m t / 


So 


1X2 


oh 64C Uj 
p 


E56C 

p 


(3.4.5) 


For turbulent f 1 ow,a=f=f ,b=0 

V 

So „ _,«?HPA*0 ,n* oftHPd' , 

'.cTfT ’ -‘MC-lTp- 1 

p V P V 


4/8 


Q 


(3.4.6) 


Q is substituted in (3.4.3) to find R. 

c h 

To evaluate the heat transfer coef f i cient , the 
correlations used for forced convection are taken as there are no 
such correlations available for natural convection. 

3.5 Formulation for dif f ernt heating cases ; 

In this section, the formulation is done for the following 
three types of heating arrangements. 

(DUniform heat flux 
(S)Triangular heat flux 
(3)Sinosidual heat flux 

The flux distributions have been shown in Figure-3. 2. For all the 
cases , eq . (3.2.3) and (3.2.4) hold ,the only change being ir \s 
source term of energy eq. (3.2.4). 

3.5.1 Uniform heat flux case : 

The governing equations are (3.2.3)and (3.2.4) which after 

being non— dimensional ised as shown before reduces to the set 
(3.4.1) with the non-dimensional parameters defined by (3.4.2). 

3.5.2 T nfetnaular flux case : 

The governing equations are (3.2.3) and (3.2.4) except for 

the source portion, energy e d * <3.2.4) being 

„ .97, G <9T,_16 /u/ , 

P C (Tr+r sq , 0<s<W/4 

w p <9 1 H &S QW,, wkxx 

=- ~q ( s-W/2) ,W/4<s<W/2 

dW max 

V/2 

In this case , power=P=tf o qnd ds 

=4/W ndq C/ W/ S ds-/ W/2 (s-M/2)ds3 

mo.* O W^- 4 r 

=W/4 ndq 

Now making non— dimensional isation as before except for 
T-Tw 

T , we get. 


1/F “=G C/ 4+B T ds-/ 2 + 2B T ds3-Q 2 

dt B it2B 
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Figure 3-z 

^ a ) Uniform heat -flux (b) Triangular hecxt -flux C^Slnosldaal heat flux 






(3. 5. 2.1) 


Sll+B") £t +Q fis 2Ms/B rO<s<B/2 

= -2M(s-B)/B ,B/2<s<B 

=0,B<s< 1+B 

=-MT r 1+B<s<1+2B 
=0, 1+2B< s<2+2B 

with continuity of T along the loop and initial condition of 6=0 
and T=0 at t=0. 

with the non-dimensional parameters given by 


F=Friction parameter=RA 1 2 /2 


G=Modified Grashof No 

M=Geometry parameters 


-S^H_ 
RhQ 2 max 

oK 

2hHA 


P C dG u 

w p oh 


_4g/9P_ 
Rhn dBG 2 


(3. 5. 2. 2. ) 


3.5.3 Sinosidual flux case : 

The governing equations are (3.2.3) and (3.2.4) except for 
the source portion, energy eq. there being 

p wVll + ! fl )=4/d q max sin(2sn/W),0<s<W/2 


In this case , powe r=P=i/ V/ 2 qn d ds=/ V ^ 2 q sin(2sn/W)frd ds 

O O max 

=dWq 

mOLK 

Now making non— dimensional isation as above except for 
* t-Tw 

T = q 7Tf ' we Set, 

man 

4 /r - dQ ~ .4tB _ __ 2 4- 2B „ . _ _2 

1/F -t=G C / T ds-/ T dsD-G 

dt » 1+2B 


1 ST Q ST 

2(1^7 St + A 5r M sin(«s/B),0<s<B 

=0,B<s< 1+B 

= -MT, 1+B<s<1+2B 
=0 , 1 +2B< s< 2+2B 

with the non-dimensional parameters given by 

2 

F=RA /2=Fr i ctional parameter 

G=Modif ied Grashof No . = — 


(3. 5. 3.1) 


(3. 5. 3. 2) 


RhdBG' 


eh 


M=Geometry parameters 


2hHA 


p C dQ L 

v p eh 


3.6 Steady state solutions^ 

The steady state equations are obtained for each case by 
making the time derivative zero in the equations (3.4.2) , (3.5.2. 1 > 


2 ? 



and (3.5.3. 1 > .There exists analytical solutions for these cases 
which are given in detail in Appendix-2 for uniform 
heating, Appendix-3 for triangular heating and Appendix-4 for 
sinosidual heating case. The results are summarised below. 

For uniform heat ing , eg . set (3.6.1) are given as 


Q=1 

T(s)=Ms+MB 

MB__ 

~1-exp( 


exp (-MB ) 

1-exp (-MB) 


-MB) 


,0<s<B 
, B<s< 1+B 


"l-exp(-MB) 


1+B<s< 1+2B 


(3.6.1) 


— , 1+2B<s<2+2B 

l-exp(-MB) 

For triangular heating, eq. set(3.6.2) are given by 
G= 1 


T(s) 


_Ms 


_ +K . , 0<s<B/2 

fc> I 

= -2M/B (s/H-B)s+K 2 ,B/2<s<B 
=K 4 ,B<s<1+B 


=K 3 exp(-Ms) , 1+B<s< 1+2B 

=K ,1+2B<s<2+2B 

. u _ DM/ „ e xp (-MB) 

where K 1 -BM/2 i_ exp( _ M §) 


K =BM/2 1 - 

K 3 l-exp(-MB) 

Kg=K MB/2 

K.=K .+MB/2 
4 1 

For sinosidual heating, eq. set 


(3.6.3) are given by 


(5=1 

T ( 5 ) =K „ — MB/n cos(ns/B) ,0<s<B 
1 

=K„,B<s< 1+B 
2 

=K 3 ex p ( -Ms ) , 1+B< s< 1+2B 

=K , 1+2B<s<2+2B 
.. . expCM(1+B)3 

where K 3 =2MB/n i_ixp7=MB) 

K 4 =K 3 expC-M(1+2B)3 
K 2 =K 3 expC-M( 1+B) 3 
K^=K .+MB/n 

*1 4 • « 

3.7 Finite difference formulation f or transient analysis 1 

For solving the transient equations, we use finite difference 

method. Implicit scheme has been used for the time derivative part 

and upwind scheme, here , backward difference, is used for the space 



derivative part. The integro-dif f erential equation is solved by 
using Simpson's rule for the integrals .The finite difference 
equations for the three cases are given below. The discretisation 
of the loop is shown in Figure- 3 . 3 . 

3.7.1 Uniform heating case: 

Let k r land n represent the no. of nodes in horizontal 
di rection,no. of nodes in vertical direction and total no. of 
nodes respectively. Source is nodalised from 1 to k ;hot segment 
from node (k+ 1 ) to (k+ 1 - 2 ) ;sink from node (k+ 1 - 1 ) to ( 2 k+l- 2 ) and 
cold segment from ( 2 k+l- 1 ) to n=( 2 k+ 21 -i ; . 

Now , taking eq. (3.4.1) and discretising as said above, we get. 


t+4 ~t k+l-4 

1 /F -77 =- 3 - a. T. - E 

v=k i.*2k+l-2 


At 


a. T t + *-a T t + *3-Q ?+ ‘ 

t L i 1 


t * _t-kl 

1 *■ * i -4 M „ , 

2~(1+B‘) At Q As Mr 1-2 to k 


= 0 ,i=k +1 to k+ 1-1 
=-MT t+ * , i=k+l to 2k+l-2 

X. 

= 0 , i=Hk+l -1 "to n 


-pt+i — .t +4 

1 4 _t + 4 4 1 r, _ . „ 

2( 1+B) "At Q " As “Or 1-1 


FGAsAt r k+l ^ + * + * + 

Q + — = L E a. T. - E a. T. -a T 3 

3 “ 4 l L \. 4 4 

Thus G l + ‘ ^ i-awwa ( 3 . 7 . 1 . 1 ) 


1+FAtQ 


t +4 


tr4 rj .2(1+B)Atr* ^ 2 ( 1 +B)At&* --“tfl 




T C1 + 

t As 


As 


L - 4 


=T t +2( 1+B)MAt,i=2 to k (3. 7.1. 2. b) 

t 


=T l , i=k+1 to k+1-1 (3. 7.1. 2. c) 

i 

, i=2k+l-1 to n (3. 7.1. 2. e) 


,2<1+B)At 

T C 1 + rr 

x. AS 


JL+4 


+2(1+B)M At3 


+ i2( 1+B )A t^t+4 

t l _ 4 is Q 


= 1 ^ , i=k +1 to 2 k + 1-2 
i 


-t + * 1+B)At n t+4 

T, C1 + — 77 « 

l AS 


2( 1+B)At rt t+4 t+A 

Q x 

As n 


(3.7. 1. 2. d) 
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=T t +2( 1+B ) MA t , i=1 (3. 7. 1.2. a) 

t 


where a =1 for i=1 , k , k+1-1 ,2k+l-2 

V. 

=4 for odd k-i and 2k+l-2-i 
=2 for even k-i and 2k+l-2-i 
3.7.2 Triangular heating cases 

The finite difference equations are exactly same as in 
(3.7.1) except at the source portion .Hence in place of (3. 7. 1.1), a 
new set is obtained f rom(3.5.2. 1 ) .The finite difference form for 
source portion gives. 






2( 1+B) 


At 


-+Q 


i. v. 

AS 


2M 


B 


< i-1 )As , i=2 to (k+1)/2 


=- ~C(i-1)As-B3,i = (k+3)/2 to k 

D 

Thus eq . (3.7. 1 .2. b) is replaced by 


T tti rj . 2 ( 1+B )At„t+i 2 ( 1+B)At 
\ C1 +"Ai Q 3 As ° V* 


— — — i— — — — ( i — i ) As , i=Z to (k+1)/2 — (3.7.t.2.b1> 
i. B 

=T t ( i_i )As-BD , i=( k+3 ) /2 to k(3.7.1.2.b2) 
i B 

3.7.3 Sinosidual heating cases 

The finite difference equations are exactly, same as in 
(3.7.1) except at the source portion. So in place of (3.7. 1 -2.b) ,a 
new equation is framed up from (3.5-.3.1). 




T t+1 -T t+ ‘ 


1 'i ‘ — — -— -=M sin to k 


2( 1+B) 


At 


-+G 


As 


n ( i— 1 ) As . 


B 


Thus eq. (3.7. 1. 2. b) is replaced by 


_ti-i r , x 2 ( 1+B )At.tti .Sil+BlAt t +1 
\ C 1+— XI Q 1 As U 'i--i 

=T l +2( 1+B)MAt sin ------- ,i=2 to k(3.7.3.2.b) 

i D 

3.8 Solution of finite difference ecij^si 

The finite difference equation set represent a coupled 
non-linear equation set in Q t+ ‘and T t+ ‘.To solve them, first 
energy equation is solved for knowing T^.Then the value of 
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\ are substituted in momentum equation to find d +l where Q fc is 
known . 

The energy equation set except for i=1 represent a bidiagonal 
matrix with subdiagonal and diagonal elements . Initial ly to make 
the system linear, Q t + ‘=Q t and T^ 4 =T^ .Knowing Q t **,T t+1 we find 
T. for i = 2, 3, • « .n . Af ter knowing T t + i ,T l + *is found from the 

1 r» 4 

energy equation with i=1. 

Now after getting all the T s at (t+1),we solve the momentum 
equation for Q^.The guessed values of Q t+4 and T^are compared 
with the fresh values for convergence within the required 
tolerance .The process is repeated till we get the converged 
solution for (t+Dth time step. 

Marching in the time direction is done till the steady state 
values are obtained, i . e . , the error between the values of G and T 
at two consecutive time steps agree within prescribed tolerance. 
3.9 Results and Discussions 

The results are shown in Graphs 3.1 to 3.18 as the plots of 
flow rate vs. time ; temperature distribution at different times 
in the loop j source temperature, hot segment tempe rature , s ink 
temperature and cold segment temperature vs. time for each of the 
heating cases. The flow rate curves show almost monotonous increase 
till steady state. The temperature vs. time curves show some 
initial oscillations before finally coming to a steady 
state . Initial ly, when there is no flow, introduction of uniform 
source results in uniform rise in temperature in the source 
portion and triangular and sinosidual sources result in 
temperature being higher at middle and lower at ends for source 
portion. For this initial per iod, conduction is the mode of heat 
transfer.As time advances, due to buoyancy, flow is 

established. Now, relatively cold fluid enters the source and hence 
the temperature in the source portion decreases. To explain the 
attainment of steady state, let us consider a hot packet of fluid 
emerging from the source and a cold packet of fluid from the sink. 
At initial time, when flow rate is low, the hot packet passing 
through the sink suffers larger temperature loss whereas the cold 
packet passing through source gains larger temperature.Now when 
the hot packet again comes to enter the source and the cold packet 
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to the sink, flow rate has increased. This results in less 
t erripe rat ur e gain for the hot packet and less temperature loss for 
the cold packet. The cycle repeats till steady state is achieved.As 
can be seen from the results, the initial non— uniformity in the 
source portion in triangular and sinosidual source cases gradually 
vanishes as natural convection is establ ished.The temperature of 
the fluid increases from entry to exit in source portion. Similar 
explanation can be given for the other portions. 
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CHAPTER-4. 

MODELLING WITH HEAT EXCHANGER 

4.1: Introduction 

In the previous chapter, we have considered the rectangular 
thermosyphon with a constant temperature heat sink for simple 

• i , *- 

analysis. In reality, the heat of the working fluid is transferred 
to a coolant in a heat exchanger .The heat exchanger may be a 
parallel flow heat exchanger where the working fluid and the 
coolant flow in the same direction or a counter flow one where the 
working fluid and the coolant flow in the opposite direct ion .The 
counter flow heat exchanger is an effective one for heat 
transfer .Here we consider the same case as the previous chapter 
with a counter flow heat exchanger as a heat sink. The analysis 
here has been done for the uniform heat flux case. However , this 
can be extended to the triangular and sinosidual heat fluxes as 
well. The thermosyphon loop for analysis has been shown in 
Figure-4. 1 . 

4.2:Formulat ion 

Here as in the case of Chapter-3, continuity 
eq. (3.2.2) , momentum eq. (3.2.3) and energy eq. (3.2.4) hold along 
with the boundary conditions (3. 2. 5). The only change here being 
T .temperature of coolant, in place of T , constant wall temperature 

C * 

in the sink portion of energy eq. (3.2.4) and in the initial 
condition (3.2.5) .The coolant temperature T^ in the heat exchanger 
is not constant.lt varies from the inlet to the outlet. Thus we 
have essentially two eq.s (3.2.3) and(3.2»4) with three variables 
Q,T and T .Thus we need another equation or to express T in terms 

e ’ 

of other two variables T,G and some known parameters such as 

coolant inlet temperature T^ .coolant flow rate G^ . 

i 

To obtain the above reduction, we consider the sink 
portion, i.e.,(W+H)/2<s<W+H/2. Considering a section from s and 
(s+ds) in this portion, we get, 
heat lost by working f luid=— pQC^dT 
heat gained by coolant=-m e C pe dT^ 
heat t ransf erred=hnd(T— T )ds 

e 

Equating all these quantities ( 
dQ=-pQC dt=-m C dT =hnd(T-T > ds 

p c pc c c 
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Now,d(T-f — •+ • ) 

p(*(. in C 

p t pc 

d(T-Te) -1_ 1 

” T- Te hnd( pQCp S~C >ds 

c pc 

To find the constant C,we have the boundary condition at 


s=W+H/2, T =T . T=T 

c e «■ «=v+H^a . 

Thus , C= <T. . HX2 -T^ . ) ex pC-tm d ( W+H/2 ) < -= Sr >^C~ ' 

p * pc 


■)3 


Substituting this value of C,we get, 
T-T = ( T 


„ „ T . )expChfid(--~ — + — - — ) ( s-W-H/2 ) 3- 

a = V + H/2 CL r pQ C m C 

p c pc 


•(4.2.1) 


Substitution of this equation in energy eq. (3.2.4) for sink 
portion gives 


^ r dT, 4h. T _ . .. 

P C ( aT + 7\ aZ ' ~~ — j(T -T . )expChrrd(- r 

w p St A Ss d »=vtH /2 ci K pQC m t 

p c pc 


-1 — + — ~ — ) ( s-W-H/2 ) 3 


Now, we have eliminated the variable T by known parameters T , 


m 


C and unknown parameters Q and T.Thus we have two 

pc 

(3.2.3) and(3.2.4) with two variables Q and T. 


equations 


4. 3? Non-dimensional i sat ion 

The non-dimensionalisation is done as in section 3.4 of 

Chapter-3.The only difference here is the non-dimensionalisation 

of T as T =-”~.Thus the non-dimensional ised equations are 

(3.4.1) except for the sink for which energy equation is 

577157 5T +q ?-=“MT expChndH/2 (-“-r=+ rri — )(s-1-2B)H 

2(1+B) St S s a=if2B r pC QQ . n» C 

— eh e pc 


(3.4.1) 


by 


1+B<s< 1+2B 

4. 3s Steady state solution 

The steady state equations are obtained for each case 
making the time derivative zero in the eq . (3.4. 1 > -There exists 
analytical solution for this case which are given below. 

The steady state equations derived from (3.4.1) are 
G LT 1 + B T ds— / 8 + ZB T ds3=G 2 (4.4.1) 

» 4 + 2B 


a dT -m 

a 3i =m 
=0 


,0<s<B 
,B<s< 1+B 


=-MT 


ft a ± + 21 


expChndH/E(- * r VvV » Cs-I-HB^^I+EB 

p eh c pc 
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=0 , 1+2B<s<2+2B 

witho continuity of T along the loop. 

This coupled equation set is solved by first solving (4.4.2) 
with Q as parameter .After getting the steady temperature 
prof i le , (4.4. 1 ) is solved for Q and then substituted in the 
solution of (4.4.2). 

Solution of (4.4.2) is given as 



exp ChndH/2(- — ) (s-1-2B)3» 1+B<s< 1+2B 

p C QG . m L 
p eh c pc 


=K 4 ,1+2B<s<2+2B (4.4.3) 

The constants are found from continuity condition of T(s). 

At s=B, K =MB/Q +K . 

2 _ 1 1 1 

s=1+B,K z =K 3 -M/a T > = ^ aii i;hndH/2(- pcjaaj' TIT’ 3 

ex pC-hTT dHB/2 ( - ~5o“~ + iTf" ’ ^ < 4 - 4 - 4) 

p ch c pc 

r — 4 -1.C7K U =U —M/n T ThrrHI-l/Pf— rr + 7C >J 


s=^ZB,K 4 =K 3 -M/Q T.^^.ChndH/Et- pra5'“ + ^c; 


s=Z+2B,K =K. 

4 1 

Now substitution of(4.4.3) in (4.4.1) gives 
G(K„-K. )=GMB/Q CFrom (4.4.4) 3=G? 

mm jt/g 

So,G=(GMB) 

Making substitutions (4.4.2) and (3.3. 10) r Q=1 (4.4.5) 

Solution of (4.4.4) set for ,K 4 and substitution of 

these constants in (4.5.3) give the steady state temperature 
distribution. 

4.5sFinite difference formulation for transient analysis 

For solving the transient equations, we use finite difference 
method.Expl icit scheme has been used for the time derivative part 
and upwind scheme, here backward difference, is used for the space 
derivative part. The integro-dif f erential equation is solved by 

using Simpson’s rule for the integrals. 

The finite difference equations are given below. 
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1/F 


Q t+1 -l? QS.s k ; l -‘ t t t t l 

" _ At —3" c z a i T i _ £ a i T i" a 4 T i :, - Q 

i-=k i*2k+t-2 14 


T 1+1 -T i TST 1 


2ll+B) At, +Q i-2 to k 


L hndH/2 , c w oc 

where z=--r and r=-^- E -- — 

m C oQ c 

e pc ^ eh 0 


m C 


=0, i=k+1 to k+1-1 

=— MT ak + ^ _ a expCz < 1-r/Q 1 ) { ( i-1 )As-1-2B>3 

i=k+l to 2k+l— 2 

=0,i=2k+l— 1 to n 


T A+4 -T t T 4 -! 4 


21 1 +B") At +Q 'As" " 1 

Thus,Q t + ‘=Q l (1— FAtQ l )+-^|--- C £ aT 1 - \ 

i=k i.=2k+l-2 


a T.^-a T t D-(4.5.1) 

V L 11 


t‘*‘=Z < 1+B >AIM+t‘ C 1igi|_ l At. T t ^ i=2 t0 k(4 . 5 . E . b) 

= T ;ti- §<Jip l -At ;i+ g<i||.d.At T; _ tii=k+1 to k+l _ 1(4 . 5 . 2 . c) 

, i=2k+l-1 to n(4.5.2.f> 


=T t [;1 _ 2( 1 +B)AtQj ) 2 ( 1+B)At(i 


a - ■ A .. -2C1+B)MAtT , , 

As AS L-l 2k+l~2 

expCz(1-r/Q l >C< i-1 )As-1-2B>3 ,i=k+l to 2k+l-3 (4.5.2. d) 


=T 4 C1- -------- — 2 ( 1+B ) AtMexpCzd-r/Q 4 ) C ( i-1 )As— 1— 2B>D3 


+ ?11±§2AtQ_ T t , i=2k+l -2 ( 4 . 5 . 2 . e ) 

As l-i 


t t 

T t+1 =T l Cl- , i=1- 


(4. 5. 2. a) 


t t~ As As 

As explicit scheme is suspectible to instability ,we choose 

the time interval at each step so as to satisfy the stability 

criterion given below. 

1-FAtG 1 > 0 
„ 2( 1+B)At 


As 


-- G >0 


1- Q -2C 1+B)MAt > 0 

As 


3 £ 
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Graph 4-i 



SOURCE TEMPERATURE VS. TIME 



□ 


Graph 4a 
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COLD SEGMENT TEMPERATURE 



a 


Graph 






4 . 6 1 So 1 ut i on of. Finite Pi f f e r en ce equat i ons 

The finite difference equations here are easy to solve as we 
have used explicit scheme. Here actually we do not need solution of 
any equations. Merely by knowing the values of Q and T at any time 
1 , we find the values at <t+1)-th time step from 
eq . s ( 4 . 5 . 1 )— (4. 5 . 2 . f ) . The time steps are chosen from the stability 
criterion at each step. 

The marching in time is done till steady state values are 
obtained, i . e ., the error between the values of G at two consecutive 
time steps agree within the prescribed tolerance. Provisions are 
also made to examine oscillatory solutions if there exists no 
steady state. This is done by allowing the marching for some 
maximum time. 

4-7sResults and discussions 

The results have been shown in the plots Graph 4.1 to 4.5. 

The graphs of flow rate vs. time ; sour ce , hot segment , s ink , cold 
segment temperatures vs. time are shown. For each portion the 
points at entry, exit and middle has been taken into account. The 
results clearly indicate oscillations. The flow rate and the 
temperatures oscillate indefinitely without finally coming into 
any steady state values. 

For each portion it indicates that the oscillation is 
convected along the portion considered. The oscillation initiates 
from the entry point and gradually passes over to the exit as 
shown from the time lag between the curves at different 
sections. The time lag between the flow rate and the temperature 
may be the cause of the observed instability. 
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CHAPTER-5 

TWO DIMENSIONAL ANALYSIS OF THE THERMQSYPHON PROBLEM 
5. ±i Introduction 

In the previous chapters, we have analysed the rectangular 
thermosyphon problem in one dimension with the coordinate running 
along the loop. However T this one— dimensional analysis leads to 
certain drawbacks. 

(1) The actual flow has been observed to be 
three— dimensional .The one— dimensional analysis gives rise to 
discrepancies between experimental and theoretical results upto 
+30X L63. 

(2) In one-dimensional analysis ,some correlations have been 
used for the friction factor f and the heat transfer coefficient 

h. The correlations are valid for forced convection 
case .Howe ver r due to their non-availability in natural circulation 
flows, we have used the same correlations. 

To avoid these difficulties, we attempt a two-dimensional 
analysis of the problem in this chapter. This eliminates the 
necessity to use any correlat ions . 

5.2 Formulation s 

For the two-dimensional formulation of the problem, the 
following assumptions are made. 

(1)The fluid is incompressible. 

<2)The flow is two-dimensional , the coordinate s running 
around the loop and r being in the radial di rection .The & 
direction is eliminated due to axi-symmetry . 

<3)The fluid properties are assumed to be constants. 

(4)Boussinesq approximation is val id, i . e ., density is assumed 
to be constant in all the conservation equations except for the 
body force term in the momentum equation where linear variation of 
density with temperature is assumed. 

i. e,p= p Cl- ft ( T-T 

a o 

(5>The velocity in the radial direction is small compared to 
the velocity in the axial direction and hence neglected. 

(6) Axial conduction and viscous dissipation are negligible 

compared to the convected energy. 

The governing equations in cylindrical coordinate are* 
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Continuity: 


dv _ 

* s -0 T i • e • w v=v ( r T t ) 
i - e Velocity is a function of ti 


me 


and 


radial 


direction. 


d v 


Momentum = p = — pg cos# 


do ,d v 1 dv. 

“"•Hi ( h— — — ) 

ds j. 2 r dr 
dr 


Integration along the loop eliminates the pressure term and 

,2 


gives 

ii = * dz -* L <~H I?' 

dr 


(5.2.2) 


Using the relation <5. 1.1), we get, 

p T 1^(~+ 3 S -> 

*odt L ^ .2 r d r 

dr 

Energy: P 0 C p (|| + v|I)=k |I) <5.2.3) 

dr 


Here axial conduction and viscous dissipation have been neglected. 

The equations (5.2.2) and (5.2.3) involve two variables T and 
G.They represent a coupled partial integro-dif f erential equation 
set. Velocity v is a function of r and t whereas temperature T is a 
function of r,s and t.So, we need initial conditions and boundary 
conditions. 

l.C.:At time t=0, T=T ,v=0 

o 

B.C . : i )Cont inuity of temperature along the loop 


ii)At r=a , v =0 
i i i ) At r = 0,^=0 


dT_ 
dr ^ 

, 0 <s<W/ 2 ( Source ) 


o 

II 

HI ^ 
<&!<$> 

,W/2<s<(W+H)/2(Hot 

segment ) 

=0 

, W+H/2< s< W+H (Cold 

segment ) 

T=T 

, ( W+H )/2<s<W+H/2( Sink) 


v ) At r= 0,~=0 

5. 3 : Non-dimens ionali sat ion 

The above equation sets are non— dimensionalised as* 

t*=tv /L s*=2s/H r*=r/a 

ch * T _ To * 

T =- — 7 — v =v/v 

qa/k 

v “Characteristic velocity=Steady state velocity for 
eh 


ch 


where 

one— dimensional thermosyphon case with constant heat flux 
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S V h 

*v* 

ST V eh qa sT 

St L 

dt* 

kL at . 

S V_ V e h 

j» * 

£V 

ST q ST* 

Sr a 

Jb * 

Sr 

tr k ,r* 

if v eh 

si* 

i f q a i 2 * 

Sr a 

Sr * 1 

B r 4 ka b r * 2 

ds=H/2 

ds* 

ST = Eqa dT^ 


Ss kH *s 

Substituting these relations in<5.2.2)— (5.2.4) and dropping 
the * symbol, we get the following non-dimensional set of 
equations . 

Momentum: ~=Gr L S ‘ + ® T ds-/®**® T dsD+IRe (---+^ (5.3.1) 

S~ r 


Energy: |I+2( 1+B) v |I=IRa(^-W |I) 

Sr 

with the boundary and initial conditions 
I .C . : At t=0, v=0,T=0 

B.C . : i )Cont inuity of T along the loop 
d v 

i i ) At r =0,~=0 


(5.3.2) 


i i i ) At r=0,~=0 
cf r 

i v ) At r=1,v=0 

v ) At r=1,f-=1,0<s<B(Source) 

or 

=0,B<s< 1+B(Hot segment) 

=0, 1+2B<s<2+2B(Cold segment) 
T=0, 1+B<s< 1+2B(Sink ) 

with the non-dimensional parameters defined by 

Gr=Grashof No — 2 — 

t.r brasnot No. 2k ^ 

eh eh 


(5.3.3) 


IRe=Inverse Reynold No.= 
IRa=Inverse Rayleigh No. 


-at . a . 4> 


P a v , 
r O eh 


p C a v . 

*0 p eh 


5. 4: Steady state analysis 

The steady state equations are derived from (5.3.3) by making 
the time derivatives zero. Thus, we get, , 2 

MomentumsGr C/* + ® T ds-/***® T dsD+IRe (^~r + ~ — (5.4.1) 
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(5.4.2) 


2 

Energy s2( 1+B) V"=IRa (-—-+— 

os *2 r or 

Or 

with the boundary condition 

i) Continuity of T along the loop 

i i ) At r=0, f-=0 

o r 


i i i ) At r=0 , ^-=0 

r 

i v ) At r=1 , v=0 (5.4.3) 

v ) At r=1 , ”=1,0<s<B 

=0,B<s<1+B 
=0, 1+2B<s<2+2B 
T=0, 1+B<s< 1+2B 

5.5iSolution for steady state 

The compexity involved in the steady state equations 
eliminates any chance of an analytical solution. The only way to 
solve these equations is by numerical means. Here we use the finite 
difference method to solve the equations .For the finite difference 
f ormulat ion , we use, central difference for the second derivatives 
and backward difference for the first order derivatives. For the 
integration encountered in the momentum equation, the Trapezoidal 
rule has been used as the intervals are of non-uniform size due to 
the presence of bends at the corners. 

The discretisation used has been shown in Figure— 5. 1 .Due to 
symmetry, without taking the total cross-section , only half of it 
has been taken for analysis. Here l,k,n and m denote the number of 
nodes in horizontal direction, in vertical direction, total number 
of axial nodes and number of nodes in radial direction 
respectively. Thus the source is nodalised from 2 to Ijhot segment 
from 1+1 to l+k-2;sink from l+k-1 to 21+k-2 and cold segment from 
21+k-1 to n=21+2k-4. 

The finite difference equations are given below. 

For interior nodes, i .e. ,j=2 to m-1 


l+k-t 

' Gr/2 C 2 a/ : T. 
i*l 


n 

E 


L J V ' J t=2l+k-2 
+ 


V . - 2 v. +v. 

a. T. . -a T H+IRe 

X. , J X. , ^ i # » J * * 


A r 


-i V - V a 

1 _1 — l-i3=0 
r i Ar 


4i 




-T. 


-2T. . +T. 


T. -T. . 


2 ( i+B ) v , L_i^i =IRa c _ix_i- 1 Iii.i_Iix.i- 1 


As. 


M 


Ar * r i Ar 

i=2 to n 

-2T. ,+T. . , T. -T. 


3, 


'•'j Ar ' j Ar 


, i=l and 21+k-2 
, i=l+k- 1 and 1 


where a. .=As. 

=As. . 

=As. ,+As. . ,for other i*s 

WJ 

=As. ,+As . F i=21+2k-4 

and As =Bendlength/2-wz/4( j-1 )Ard/H for i=1,2,l,l+1, 

l+k-1 , l+k,21+k-2 and 21+k-1 

=As for other i*s 

For j=1 , i . e . ,at r=0,“=0 and f-=0 

Or Or 

"t ^ 1 OT 2 ^ 

The terms — ” and - -z- become — — and — -by L'Hospital *s rule, 
r or r Or A 2 ^ 2 

Sr Sr 


These conditions also give v =v and T. =T 

1 o v , 1 i,a 

During discretisation of second derivatives r we get, for 
j=1 , expressions involving v a nd T. for which the above 

O L ,0 

substitutions are made. 

Thus, we get, for j=1 


l-rk-i „ v-v 

* Gr/2 C X a. T. a. T. -a T 3+2IRe ; =0 

1,1 t. , 1 l = 2 1 1 k-2 l , 1 L , 1 1,1 1,1 * 2 

u=l A r 

T . -T . T -T . 

and ( 1+B) v i=2 to n 

i As. . 2 

Ar 

T -T T -T 

( 1+B) V i = 1 

‘ As .,* Ar 8 


At r=1 , i .e . , j=m, 
v=0,or,v^=0 

”=1 ,0<s<B,or,T. ,=T. . +Ar,i=2 to 1 

S r <-rJ 

=0,B<s<1+B,or,T t =T^ ^_ 4 »i=l+1 to l+k-2 

T=0, 1+B<s<1+2B,or,T ^=0, i=l+k-1 to 21+k-2 

™=0, 1+2B<s<2+2B,or,T. =T. . ,i=21+k-1 to n and 1 

Sr , j v , j-4 

Rearrangement of the above equations give. 
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Energy equation set(5.5.1) is given by 


IRaAs. IRaAs. 

For j = l:T Cv + ~-^3 —j T . =T. . ,i=2 to n 

' J J (1+B)Ar (1+B)Ar + 1 


=T . i i=1 

n,j 

IRaAs. . IRaAs. . 

For j=2 to m-lsT C 1^1 ( _ 1+ 2 j 3+T Cv + --^(2 — A-r 

' J 2(1+B)Ar J 1 v ' J J 2(1+B)Ar J_1 

IRaAs. 

+T. r --i-3=T. 7 i =2 to n 

' J 2(1+B)Ar 


)□ 


=T .,i=1 

n,J 

For j=m=T ,=T. +Ar,i=2 to 1 

L , j t , J - 1 

=T. . ,i=l +1 to l+k-2 

L , j-i 

= 0 , i=l+k -1 to 21 +k -2 

=T. ,i=21+k— 1 to n and 1 

L , J - 1 

The momentum equation set(5.5.2) is given by 


n- - . . GrA r 

For j = 1 : v . -v . =^ 75 - 

j j+i 4IRe 


2 btk-A 


C E a. T. ? , a. T. -a T 

l,J l.J 1 . =2 l + k-2 L , J i,J 1,0 


j-Z to m* v . C-l+r—H+v . C*----3-v . 

j-i j-i j J-1 j+i 


GrAr 


2 l+k-l 


[ 2 a T - ,1, a. T. -a T 

2IRe i,j V,j i = 2l+k-2 , J L,J l,j l,j 


i=l 


j =m s v . =0 

J 

5.6sSolution Techniq ue 

•Jhe above equation sets (5.5.1) and (5.5.2) represent a 
non— linear set of equations. Hence they are to be solved 
iteratively. Each set represents a system of tridiagonal matrix 
with sub-diagonal * diagonal and super-diagonal elements. 

To solve the set of equations , ini tially some guess values of 
the vectors fv>° and IT>° are assumed. Fi rst , the tridiagonal set 
(5.5.1) is solved. The quantities v. and T. appearing in the 
elements of the sub-diagonal , diagonal , super-diagonal or RHS 

vectors are taken to be v. and T^ . res pectively .Thus we get the 

1- ^ # J 

value of CT } 4 for first iteration step. Then the tridiagonal 

set ( 5 . 5 . 2 ) is solved. The quantities v ,T appearing in the 

J ** * J 

elements of the sub-diagonal , diagonal , super-diagonal or RHS 
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vectors are taken to be v°and T . respect ively. Thus we get the set 

i . 4 J 4 L ' 3 

{TJ .Now the sets -Cv} and tT> are compared with the previous set 
values for convergence. If the system has not converged, then the 
new values are taken as guess values and the process is repeated 
till convergence. Provisions are also made to detect if the 
solution ,in stead of converging , osci 1 lates .For this purpose, 
number of oscillations is allowed upto a maximum and the values 
are noted at each iteration to check if oscillations are present. 
5.7sResults and discussions 

The results are plotted in Graphs 5.1 to 5.5 
as velocity vs. iteration numbersand temperature along the source 
,hot segment, sink and cold segment vs. iteration numbers. The 
curves show the oscillations in the values of the above quantities 
with iterations .From the oscillatory behaviour, it can be 
concluded that there exists no steady state for this case with the 
chosen values of the parameters .Any small deviation from the 
steady state results in oscillations. This agrees with the observed 
instabilities or even flow reversals in the rectangular 
thermosyphon loop. 
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VELOCITY VS. ITERATION 



D 


ITERATION WJM3ER ^ 

v«l(l) + v«l(2) ❖ vel(«3) A vel(4) x v«l(5) 



SOURCE TEMPERATURE VS. 




HOT SEGMENT TEMPERATURE VS. 



2aniva3dH2i 






CONCLUSION 


The analysis of the internal natural ci rculat ion flow in the 
rectangular thermosyphon loop has been done starting from the 

simple case of a thermosyphon loop with point source and point 

sink to the two-dimensional analysis of the loop with extended 

source and sink with the heat exchanger as the sink. It has been 

assumed throughout that there is no initial flow in the loop. The 

steady state, transient and stability characteristics of each case 

has been analysed in detail. 

It has been shown that the behaviour of the loop depends 
largely on the prescribed model. The one-dimensional rectangular 
thermosyphon loop with constant temperature heat sink shows no 
instability, whereas the same loop gives rise to instability with 
the introduction of the heat exchanger .The two dimensional 
analysis also leads to instability. 

The loop configuration is one of the few alterable factors in 
the thermosyphon design. The geometrical dimension of the loop can 
be changed so as to get rid of the operational and safety problems 
caused by the periodic instability in the natural convection 
loops . 

However, the actual process occuring in the nuclear power 
reactors is much more complex than the case analysed.lt comprises 
of a number of such parallel loops with the pressure drop varying 
largely at different sections .Again , after the failure of the 
elecrical power supply, due to the high inertia of the flywheel of 
the primary pump, flow is maintained for some period of time and 
thereafter natural circulation takes over. This transition from 
forced flow to natural circulation flow has not been taken into 
account in our case. The improvements in the modelling so as to 
include the parallel loops, consideration of the pressure drops at 
different sections of the loop and the initial condition of 
transition from forced flow to natural circulation flow will throw 
more light in the design of the primary heat transport system of 
the nuclear reactor so as to remove the decay heat effectively to 
ensure safety of the plant. 
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APPENDI X-I 


Let us consider a fluid particle passing through the source 
at temperature AT for a time At. Let T. and T represent the 

temperature of the fluid particle before and after passing through 
the source and Q the flow rate. 

Now, energy equation gives 
k<AT— T> 

are considering a particular particle 
mot ion , space , i . e ., posit ion of the particle becomes a function of 
time and hence in all we have one independent coordinate t. 

So, we get,— +kT=k AT 
Now, T =A exp(-kt) 

OX 


SJ 

St 


G *T 
A S~s 
As we 


T =AT 

p»- 

So, T(t)=A exp(-kt)+AT 

The boundary conditions are at t-0 , T=T^ 

at t=At , T=T 

out 


So , T =A+AT 

tri 

T =A exp<— kAt)+AT 
°’ ,Jt LA| -kAt 

and T -T. =A(-1+e ) = (AT-T. ) ( 1-e ) 

out i r» t r» 

Similarly, for heat sink, 

T -T =( -AT-T. )<1-e" kAt ) 

out tr> 
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APPENDIX -I I 


STEADY STATE SOLUTION FOR UNIFORM HEAT SOURCE 
The steady state equations derived from (3.4.1) are 

... 2B"’ “ o 

GDf B T ds-X 4 + ““T dsO=Q 2 (A2.1) 


Q ds =M ' 0<S<B 

=0 , B< s< 1+B 

= -MT T 1+B<s< 1+2B (A2.2) 

=0 , 1+2B<s<2+2B 

with continuity of T along the loop. 

This coupled equation set is solved by first solving (A2.2) 
with Q as parameter .After getting the steady temperature profile, 
(A2.1) is solved for G and then substituted in the solution of 
(A2.2) . 


Solution of (A2.2) is given as 
T(s)=Ms/G+K lT 0<s<B 


=K 2 ,B<s< 1+B 


(A2.3) 


=K 3 exp(-Ms/G ) , 1+B<s< 1+2B 
=K 4 , 1+2B<s<2+2B 

The constants are found from continuity condition of T(s). 

So, at s=B,K„=MB/G+K. 

HI _ 

at s=1+B, K 2 =K 3 expC-M( 1+B)/QT 
at s=1+2B,K 4 =K 3 expC-M(1+2B)/QT 
at s=2+2B,K 4 =K 1 

Solving these four equations for K ,K„,K and K_,we find, 
u _mr ex£(— MB/G) 

k 1 -k 4 -mb/g “rf;57r S B7S) 

K =MB/Q 7-S57J5T (AH. 4) 

2 1-exp(-MB/Q> 

K -MB/G expLM(1+B)/Q3 

k 3 -mb/g 1 _ exp( _ MB/ ts ) 


Now, substitution of (AH. 3) in (AH.1) gives 
G(K 2 -K 4 >=GMB/G CFrom A2.4T=^ 

so7g=<gmb> 4 ' 3 

Making substitutions (3.4.2) and (3.3. 10) ,G=1 (A2.5) 

The temperature distribution is found by substituting (A2.4) 
and (A2.5) in (A2.3). 


4 ? 



APPENDIX -I I I 


STEADY STATE SOLUTION FOR TRIANGULAR HEAT SOURCE 


The steady state equations derived from 
GD/_ T ds-T T dsO=G (A2.1) 

b i4 as 


(3.5.2. 1) are 


G ” =2sM/B ,0<s<B/2 
as 

=-2M(s-B)/B ,B/2<s<B 
=0 ,B<s< 1+B 

= -MT , 1+B<s< 1+2B (A3. 2) 

=0 , 1+2B<s<2+2B 

with continuity of T along the loop. 

This coupled equation set is solved by first solving (A3. 2) 
with G as parameter .After getting the steady temperature profile 
(A3.1) is solved for Q and then substituted in the solution 
(A3. 2) . 


Solution of (A3. 2) is given as 
T(s)=Ms^(QB)+K 1 ,0<s<B/2 


=-2Ms/(BG>*(s/2-B)-H( 2 , B/2< s<B 


=K 4 , B<s< 1+B 


■ ( A3 . 3 > 


=K 3 exp(-Ms/G ) ,1+B<s<1+2B 
=K 5 , 1+2B<s<2+2B 

The constants are found from continuity condition of T(s). 
s=B/2 , MB/ ( 4G ) -i 
at s=B , K_+MB/G=K , 


So , a t s=B/2 , MB/ ( 4G ) +K JJ =3MB/ ( 4G ) +K 2 


'2"‘~' “ 4 

at s=1+B, K 4 =K 3 e>;pC-M( 1+B)/G3 

at s=1+2B r K 5 =K 3 expC-M(1+2B)/G3 

at s=2+2B.K 5 =K 1 


Solving these five equations for K. ,K ? ,K„,K.and K,_,we find, 
K -K -MB/(2G> 

K l“ K 5 _MB/ _ 1-exp(-MB/5> 


K 2 =-MB/(2G)+K 1 


(A3. 4) 


expCM(1+B)/GD 
K 3 -MB/(2Q) 1 _ eX p(_M B /S) 


K 4 =K 1 +MB/(2Q> 

Now, substitution of (A3. 3) in (A3.1) gives 

G ( K . -K- ) =GMB/ ( 2Q ) CFrom A3.4D=Q 2 

^ ^ 

So , G= ( GMB/2 ) 

Making substitutions (3. 5.2.2) and (3.3. 10) ,G =: 1 (A3. 5) 


of 


4 * 



and 


The temperature distribution is found by substituting (A3. 4) 
(A3. 5) in (A3. 3). 
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APPENDIX -IV 


STEADY STATE SOLUTION FOR SINQSIDUAL HEAT SOURCE 

The steady state equations derived from (3. 5. 3.1) are 
GCf‘ +I> T ds-/* + **T ds3=Q a (A4.1) 

® 1 + Sfi 

- dT 

Q =M sin(ns/B) ,0<s<B 
=0 ,B< s< 1+B 

= -MT , 1+B<s< 1+2B < A4 . 2 ) 

=0 » 1 +2B< s< 2+2B 

with continuity of T along the loop. 

This coupled equation set is solved by first solving (A4.2) 
with Q as parameter .After getting the steady temperature profile, 
(A4.1) is solved for G and then substituted in the solution of 
<A4.2> . 

Solution of (A4.2) is given as 

T(s)=-MB/(Gn)cos (ns/B)+K 1 ,0<s<B 

=K 2 ,B<s< 1+B (A4.3) 

=K 3 exp<-Ms/G ) , 1+B<s< 1+2B 
=K 4 , 1 +2B< s< 2+2B 

The constants are found from continuity condition of T(s). 

So, at s=B,K 2 =MB/(nG)+K 1 

at s=1+B, K 2 =K 3 expC-M(1+B)/G3 
at s=1+2B,K 4 =K 3 expC-!1(1+2B>/G3 
at s=2+2B,K 4 =K 1 -MB/(frG) 

Solving these four equations for K,j,K 2 ,K 3 and K 4 ,we find. 


K 2 =2MB/(nG) 


~1-exp(-MB/§> 

„„„ , . ~ . expCMC 1+B)/GI3 
K 3 =2MB/< n G 

K 4 =K 3 ex pC-M C 1+2B ) /§□ 

K .=K ,+MB/ (Gre ) 

1 4 


(A4.4) 


Now, substitution of (A4.3) in (A4.1) gives 
G (K 2 - K 4>=2GMB / (rcQ) CFrom A4.43=Q 2 
SoTG=i2GMB/n>^ a 

Making substitutions (3. 5. 3. 2) and (3.3.10),Q=1 (A4.5) 

The temperature distribution is found by substituting (A4.4) 
and (A4.5) in (A4.3). 
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